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Genetic correlations greatly increase mutational robustness and can both reduce and enhance 

evolvability 
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Mutational neighbourhoods in genotype-phenotype (GP) maps are widely believed to be more likely to share 
characteristics than expected from random chance. Such genetic correlations should strongly influence evolu¬ 
tionary dynamics. We explore and quantify these intuitions by comparing three GP maps - a model for RNA 
secondary structure, the HP model for protein tertiary structure, and the Polyomino model for protein quaternary 
structure - to a simple random null model that maintains the number of genotypes mapping to each phenotype, 
but assigns genotypes randomly. The mutational neighbourhood of a genotype in these GP maps is much more 
likely to contain genotypes mapping to the same phenotype than in the random null model. Such neutral cor¬ 
relations can be quantified by the robustness to mutations, which can be many orders of magnitude larger than 
that of the null model, and crucially, above the critical threshold for the formation of large neutral networks of 
mutationally connected genotypes which enhance the capacity for the exploration of phenotypic novelty. Thus 
neutral correlations increase evolvability. We also study non-neutral correlations’. Compared to the null model, 
i) If a particular (non-neutral) phenotype is found once in the 1-mutation neighbourhood of a genotype, then 
the chance of finding that phenotype multiple times in this neighbourhood is larger than expected; ii) If two 
genotypes are connected by a single neutral mutation, then their respective non-neutral 1-mutation neighbour¬ 
hoods are more likely to be similar; iii) If a genotype maps to a folding or self-assembling phenotype, then its 
non-neutral neighbours are less likely to be a potentially deleterious non-folding or non-assembling phenotype. 
Non-neutral correlations of type i) and ii) reduce the rate at which new phenotypes can be found by neutral 
exploration, and so may diminish evolvability, while non-neutral correlations of type iii) may instead facilitate 
evolutionary exploration and so increase evolvability. 
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I. AUTHOR SUMMARY 


Evolutionary dynamics arise from the interplay of mutations 
acting on genotypes and natural selection acting on phenotypes. 
Understanding the structure of the genotype-phenotype (GP) 
map is therefore critical for understanding evolutionary pro¬ 
cesses. We address a simple question about structure: Are the 
genotypes positively correlated? That is, will the mutational 
neighbours of a genotype be more likely to map to similar phe¬ 
notypes than expected from random chance? John Maynard 
Smith and others have argued that the intuitive answer is yes. 
Here we quantify these intuitions by comparing model GP maps 
for RNA secondary structure, protein tertiary structure, and pro¬ 
tein quaternary structure to a random GP map. We find strong 
neutral correlations: Point mutations are orders of magnitude 
more likely than expected by random chance to link genotypes 
that map to the same phenotype, which vitally increases the 
potential for evolutionary innovation by generating neutral net¬ 
works. If GP maps were uncorrelated like the random map, 
evolution may not even be possible. We also find correlations 
for non-neutral mutations: Mutational neighbourhoods are less 
diverse than expected by random chance. Such local hetero¬ 
geneity slows down the rate at which new phenotypic variation 
can be found. But non-neutral correlations also enhance evolv¬ 
ability by lowering the probability of mutating to a deleterious 
non-folding or non-assembling phenotype. 


II. INTRODUCTION 

In a classic paper m, published in 1970, John Maynard 
Smith introduced several key ideas for describing the structure 
of genotype-phenotype (GP) maps. He first outlined the con¬ 
cept of a protein space, the set of all possible sets of amino acid 
chains, and suggested that for evolution to smoothly proceed, 
these should be connected as networks of functional protein 
phenotypes that can be interconverted by (point) mutations. He 
then argued that one criterion for such networks to exist is for a 
protein X to have at least one mutationally accessible neighbour 
which is “meaningful, in the sense of being as good or better 
than X in some environment”. In other words, if X has N mu¬ 
tational neighbours, then the frequency / of “meaningful” pro¬ 
teins in its mutational neighborhood should satisfy / > 1/N. 
He pointed out that this was likely to be true in part due to the 
ubiquity of neutral mutations, which had been famously pro¬ 
posed by Kimura m and King and Jukes HI just a few years 
prior to his paper. But he also gave a second reason for ex¬ 
pecting connected networks, namely that, “There is almost cer¬ 
tainly a higher probability that a sequence will be meaningful 
if it is a neighbour of an existing functional protein than if it is 
selected at random.” This concept that mutational neighbours 
differ from the random expectation is what we will call genetic 
correlations. 

Following Maynard Smith, many authors have explored the 
role of networks of genotypes connected by single point muta¬ 
tions. Lipman and Wilbur d first showed that large networks of 
mutationally connected genotypes mapping to the same pheno¬ 
type are found in the Hydrophobic-Polar (HP) model for protein 
folding, introduced by Dill (HO. They also pointed out that 
neutral mutations allow a population to traverse these networks, 
facilitating access to a larger variety of alternate phenotypes. 
Schuster and colleagues 171 developed these themes further us¬ 
ing detailed models for the secondary structure of RNA (H. 
They coined the term “neutral network” to describe sets of mu- 
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tationally connected genotypes that map to the same secondary 
structure phenotype. As RNA secondary structure is fairly easy 
to calculate and thermodynamics based models such as the Vi¬ 
enna package are thought to provide an accurate prediction of 
real RNA secondary structure ||9j[T0]|, the nature of neutral net¬ 
works in these models has been extensively studied Cl Emu- 
da. Since these pioneering works, neutral networks have been 
considered in GP maps of other biological processes, includ¬ 
ing models for gene networks ifTTl [T^ metabolic networks m 
and the Polyomino model for self-assembling protein quater¬ 
nary structure l(20]| . 

From these studies of model systems a number of basic prin¬ 
ciples have emerged, much of which has been reviewed in im¬ 
portant books by Wagner ifTSlfT^ . Firstly, for neutral networks 
to exist, the GP map should exhibit redundancy, where multiple 
genotypes map onto the same phenotype. This many-to-one na¬ 
ture of the mappings is illustrated in Fig. [T}\. Redundancy is of 
course closely linked to the existence of neutral mutations ElEl, 
although the relationship between these concepts is not entirely 
unambiguous. In the theory of neutral evolution, a mutation 
may lead to a slightly different phenotype, but as long as the 
change in fitness is small enough not to be visible to selection, 
it is considered to be effectively neutral I2TI . Whether selection 
can act depends on the degree of phenotypic change, the envi¬ 
ronment, and other factors such as the population size and mu¬ 
tation rate. Therefore, identifying whether or not a mutation is 
neutral can be complex, and the answer may vary as parameters 
external to the GP map change with time. So while redundancy 
only couples identical phenotypes, and so is a more restrictive 
concept than neutral mutations, it has the advantage of sidestep¬ 
ping the subtle issues listed above and is therefore more easily 
applicable to the study of a static GP map. 

The second basic principle to emerge is that the number of 
genotypes per phenotype (the redundancy) can vary, leading to 
phenotype bias, as depicted in Fig.l^B. Thirdly, it is generally 
the case that the larger the redundancy, the greater the mean 
mutational robustness of genotypes mapping to that phenotype. 
Fourthly, the larger the neutral network, the greater the variety 
of alternative phenotypes within one (non-neutral) point muta¬ 
tion of the whole neutral network, leading to a positive correla¬ 
tion between robustness and measures of evolvability that count 
the number of different phenotypes that are potentially accessi¬ 
ble lEl. 

Finally, a key principle emphasised by Maynard Smith |[T|, 
but which has earlier roots in concepts such as the shifting bal¬ 
ance theory of Sewall Wright ||23l, is that neutral mutations al¬ 
low a population to access, over time, a wider variety of po¬ 
tential alternative phenotypes than would be available around a 
single genotype (uiniiisi. Evidence for the key role of these 
networks in promoting evolutionary innovation has been found, 
for example, in experiments on RNA structures 1241 [25l and 
transcription factors 12^ . 

The main focus of this paper is genetic correlations. To ex¬ 
plore and quantify how they affect concepts such as neutral 
networks, robustness and evolvability, we study genetic corre¬ 
lations in three of the GP maps mentioned above. These are 
the sequence to RNA secondary structure map and HP model 
for protein folding (tertiary structure), which have been exten¬ 
sively studied, as well as the more recently introduced Poly¬ 
omino model for self-assembling protein quaternary structure. 
Several properties of these three GP maps have recently been 
compared ||20| [271 and we summarise some of the similarities 


and differences between them in the methods section. However, 
a detailed investigation of their genetic correlations has not yet 
been considered. 

A key question we consider is how to define genetic corre¬ 
lations in a quantitative way. For that one needs some kind of 
uncorrelated null model for how genotypes are distributed over 
phenotypes with which to compare the full biophysical systems. 
We employ a model we call the random GP map. It has the same 
number of genotypes mapping to each phenotype as the bio¬ 
physical GP map to which it is being compared, as well as the 
same basic type of genotype space (alphabet size and genome 
length), and nodes (genomes) that are linked by single muta¬ 
tions if they differ by one locus. The difference is that the geno¬ 
types are randomly distributed across the genotype space. Of 
course one does not expect biological systems to have such a 
random distribution, but it is not at all straightforward to think 
of a better null expectation. The great advantage of having such 
a null model, even if one knows that it is has limitations, is that 
it allows us to quantitatively contrast how the biophysical geno¬ 
types are organised across the genetic space, which should shed 
light on nature of the correlations that Maynard Smith intro¬ 
duced. 

The paper is organised as follows. We first define our mod¬ 
els in the methods section. Next we examine neutral correla¬ 
tions, schematically illustrated in Fig.[2C, through considering 
various measures of robustness that quantify the relative like¬ 
lihood that mutationally neighbouring genotypes possess the 
same phenotype. We then perform a similar analysis compar¬ 
ing the biological GP maps to the random map for non-neutral 
correlations. Since these different kinds of correlations all mod¬ 
ulate the way that novel variation arises through random muta¬ 
tions, we finish by commenting on how correlations affect sub¬ 
tle interplay of robustness and evolvability |[l5l[28l, and also 
briefly suggest a few other forms of correlation that could be 
studied in GP maps. 


III. METHODS 

A. Biological GP maps: RNA, Polyomino and HP models 

We consider three separate GP maps for low-level self¬ 
assembling biological systems. Firstly, a model for RNA sec¬ 
ondary structure 111 that determines which bases in an RNA 
sequence form bonded pairs. Secondly, the HP lattice model 
for protein tertiary structure Eia, that determines the three- 
dimensional shape of a folded amino acid chain. And thirdly, 
the Polyomino model for protein quaternary structures 1201 [2^ 
l30l . where quaternary structure of proteins is the topological 
arrangement of separate folded amino acid chains. All three of 
these models have been previously compared in ref. ll2Qll . and 
below we briefiy outline the three different systems. 


1. Vienna package for RNA secondary structure 

In the widely studied RNA secondary structure GP map 171 
[8l[nHl6]|, the genotypes are sequences made of an alphabet 
of four different nucleotides, and phenotypes are the secondary 
structures, which describe the bonding pattern in the folded 
structure with the lowest free energy for the given sequence. 
Here, we use the popular Vienna package ISl, which uses an 
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B) Phenotype Bias 
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C) Correlations 
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FIG. 1. Schematic depiction of the GP map properties of redundancy, phenotype bias and neutral correlations. Phenotypes are represented 
by colours, genotypes as nodes and mutations as edges. A) Each colour appears multiple times with uniform redundancy. B) Some colours 
appear more often than others, demonstrating a phenotype bias. C) A rearrangement of the colours from the middle plot illustrates positive neutral 
correlations where the same colours are more likely to appear near each other than would be expected by random chance arrangement. The black 
box surrounding the six orange genotypes depicts a single component (a set of genotypes connected by neutral point mutations, also called a 
neutral network) of the orange phenotype. Such positive neutral correlations enhance the probability that such neutral networks occur. 


empirical free-energy model and dynamic programming tech¬ 
niques to efficiently find the lowest free energy structures. We 
use Version 1.8.5 with all parameters set to their default val¬ 
ues. The RNA GP map for length L is referred to as RNAL. 
We present results from RNA 12, with 4}^ ^ 16.8 x 10^ geno¬ 
types mapping to 57 folded pheontypes, RNA15 with 4}^ = 

I. 07 X 10^ genotypes mapping to 431 folded phenotypes and 
and RNA20, with 4^^ ^ 1.10 x 10^ genotypes mapping to 

II, 218 folded secondary structure phenotypes. 


2. HP lattice model for protein tertiary structure 


The HP lattice model represents proteins as a linear chain of 
either hydrophobic (H) or polar (P) amino acids on a lattice O . 
A simple interaction energy function is used between non- 
adjacent molecules in the chain. Hydrophobic-hydrophobic in¬ 
teractions are energetically favourable. We use an implementa¬ 
tion with the interaction energy E between the different poten¬ 
tial pairs being classified as E^hh = — 1, with = Epp = 0, 
as widely chosen by other authors, see e.g. 121 ED- A phe¬ 
notype is defined for each shape (fold) on the lattice that is 
the unique free energy minimum of at least one sequence. If 
a sequence has more than one structure as its minimum, then 
it is considered not to fold properly, and so is categorised as 
belonging to the (potentially deleterious) general non-folding 
phenotype. The compact HP model restricts the possible folds 
to those which are maximally compact, in an attempt to cap¬ 
ture the globular nature of in vivo proteins |[3^ . We make use 
of both compact and non-compact HP GP maps by consider¬ 
ing both the GP map for all folds of length, L = 24 (denoted 
HP24) and all compact folds on the 5 x 5 square grid (of length 
L = 25, denoted HP5x5). For the non-compact HP24 GP map 
there are 2^^ 16.8 x 10^ genotypes and 61,086 folded pheno¬ 

types, while for the compact HP5x5 model the 2^^ ^ 33.6 x 10^ 
genotypes map to a much smaller set of 549 unique folded phe¬ 
notypes. 


3. Polyomino model for protein quaternary structure 

Protein quaternary structure describes the topological ar¬ 
rangement of separate folded proteins that self-assemble into 
a well-defined cluster. The Polyomino GP map is a recently 
introduced lattice model which, in the spirit of the HP lattice 
model for tertiary structure, provides a simplified but tractable 
GP map for protein quaternary structure, as described in more 
detail in ||20|- Very briefly, it employs a set of Nt square tiles 
with Nc interface types, together with a set of rules that denote 
how the interfaces bind. These sets are specified by genotypes 
in the form of linear strings. In this work, Nc = 8 and we spec¬ 
ify that the interface types interact in ordered odd-even pairs, 
such that the following interface types interact (1 ^ 2, 3 ^ 4, 
5 6). 

The conversion of the tile set (genotype) into a bounded shape 
(phenotype) is achieved through simulating a lattice based self- 
assembly process. This begins with seeding where the first tile 
encoded in the genotype is placed on the lattice. Thereafter 
other tiles are randomly placed. If they form a bond they are 
kept, otherwise they are discarded. This assembly process is 
repeated until either no more bonds can be formed or else the 
number of tiles grows without limit. For each set of tiles, as¬ 
sembly is attempted several times. If a given set of tiles self- 
assembles into a unique bounded shape, the genotype is con¬ 
sidered to map to that shape (phenotype). If on the other hand, 
the tile set does not always assemble to the same shape, or if 
it assembles to an unbounded shape (as occurs in sickle cell 
anaemia for example), then the genotype maps to the undefined 
phenotype (UND). 

The GP map resulting from Nt kit tiles and interface types 
is denoted as In this work, we consider £' 2,8 which has 

1.7 X 10^ genotypes mapping to 13 different self-assembling 
phenotypes and the larger space which has 6.9x10^^ geno¬ 
types mapping to 147 phenotypes. All parameters used for sim¬ 
ulations were the same as in ref. 1201 . 


4. A deleterious phenotype in all three GP maps 

We also distinguish a deleterious phenotype (del) in all three 
GP maps. For the RNA GP map, this occurs when the unbonded 
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strand is the free-energy minimum, so that the strand does not 
fold. For the HP model this occurs when a sequence does not 
have a unique ground state structure, which is interpreted as 
the protein not folding. For the Polyomino GP map this oc¬ 
curs when the set of tiles produces an UND phenotype. Note 
that, depending on the environment, many other phenotypes 
may also be deleterious, but the del phenotype will always be an 
evolutionary dead end. For RNA this del phenotype makes up 
85% of RNA12, 65% of RNA15 and 33% of RNA20. In the HP 
model the fraction is typically larger, consisting of 98% of the 
HP24 map and 82% for the HP5x5 mapping, while for the Poly¬ 
omino GP map we find 54% of 82,8 and 80% of For the 
RNA this fraction decreases with increasing genotype length L 
(asymptotically the fraction fdei, of genotypes mapping to the 
deleterious phenotype scales as fdei ~ 21.4 x 0.82^ (3^ ) while 
for the Polyomino GP map fdei increases for larger system size, 
and in the HP model the trend remains ambiguous oa. 


5. Similarities and differences in the three GP maps. 

Recently, direct comparisons between the properties of these 
GP maps have begun to be considered (201 El- A common 
feature of all three GP maps is that they correspond to self- 
assembly processes of molecular structures in biological sys¬ 
tems. Thermodynamics drives the assembly process in each 
case. They differ in that the RNA and HP model have linked 
units that are fixed in size, while the polyomino model has dis¬ 
connected units that can assemble into multiple sizes. The GP 
maps share properties such as redundancy, bias and phenotypic 
robustness (20l[27l (see Fig. [^, which are seen in a much wider 
range of GP maps (TsKTbl. However, the nature of these prop¬ 
erties varies between the GP maps. For example, the polyomino 
and RNA GP maps typically have fewer phenotypes and larger 
neutral networks than phenotypes in the HP model. Another 
key variable that differs between the GP maps is the base K. 
For RNA K = 4, HP K = 2 and for Polyominoes K = S. 
These different genotypic topologies as well as the way in basic 
units interact to generate phenotypes can affect the formation of 
neutral networks. For example, in RNA, a CG bond in a stem 
motif cannot be turned into a GC bond without breaking the 
bond (351, a form of neutral reciprocal sign epistasis (361 . For 
a secondary structure stem motif of n bonds, this phenomenon 
breaks the neutral set up into around 2 ^ separate components, 
often of similar size, that are connected by point mutations in¬ 
ternally, but which need at least two mutations to be connected 
together. These different components of the full neutral set are 
separate neutral networks. With K = S and asymmetric bond¬ 
ing of bases (1^2 and 2^1 being distinct), the Polyomino 
GP map may be similarly susceptible to this form of epista¬ 
sis. The HP lattice model does not have this feature due to 
bonds being formed from non-adjacent neighbouring HH inter¬ 
actions. How these connectivity patterns scale with increasing 
L remains an open question, in part because the spaces grow 
exponentially with length, and so rapidly become much more 
difficult to comprehensively analyse. 


B. The random GP map as an uncorrelated null model 

As discussed in the Introduction, in order to quantify genetic 
correlations we must first define an uncorrelated null model to 


which the biophysical GP maps can be compared. Here we em¬ 
ploy a random GP map that was recently explicitly introduced 
for analysing whole GP map properties in ref. (141, but has also 
been used implicitly in many earlier works see e.g. da [351 [371, 
although we believe this is the first time this random model has 
been used to define correlations. 

The random GP map shares the following properties with the 
biological GP map to which it is being compared: the same al¬ 
phabet size K, genome length L, number of 1-mutation neigh¬ 
bours {K — 1)1/, number of genotypes Nq = number of 
phenotypes Np, and frequencies fp, defined as the fraction of 
all genotypes that possess phenotype p. We summarise the GP 
map nomenclature used in this paper in Table[T| which compares 
what is shared and what is different between the biological and 
the random GP maps. 

With these key global GP map properties fixed, the only dif¬ 
ference between a biological GP map and its associated random 
GP map is that the Fp = fp x Nq genotypes for each pheno¬ 
type p are each randomly assigned to the set of Nq possible 
genotypes. As phenotypes are randomly assigned, departures 
in properties between the two versions of the GP map may be 
considered to be due to correlations, that is , in contrast to the 
random GP map, the mutational neighbourhood of a genotype 
in the biological GP maps is affected, for example, by what phe¬ 
notype it maps to. We note that these correlations can be very 
complex, and depend not only on the identity of the phenotype, 
but also on higher order features such as the identity of one or 
two or more phenotypes in the direct neighbourhood (higher or¬ 
der correlations). It is almost certainly also true that, depending 
on the phenotype, the correlations may also depend on which of 
the L genomic positions is mutated (see e.g. Maynard Smith’s 
word game described in the discussion). However, in this paper 
we mainly focus on the simplest kinds of correlations; for ex¬ 
ample for neutral correlations we mainly look at effects that are 
captured by the concept of robustness. 

Other null models are conceivable. For example, in ref. na 
the authors used an approach based on network theory 
comparing the topology of neutral components to Erdos-Renyi 
networks and scale-free networks. While this type of network 
is helpful for understanding the topology of neutral networks 
themselves, a focus of this work here is not just how genotypes 
of the same phenotype connect to each other, but how pheno¬ 
types are arranged in relation to each other. Besides their sim¬ 
plicity, an advantage of using the random GP maps for this pur¬ 
pose is that the overall connectivity of the genotype space is left 
intact, along with several global properties of the map, allow¬ 
ing the way phenotypes are arranged to be directly considered. 
Moreover, this random map has been used (implicitly and ex¬ 
plicitly) throughout the literature, see e.g. CSEIIISI, and so it 
is of general interest to carefully analyse some of its properties. 


IV. RESULTS 

A. Phenotypic robustness and neutral correlations 

The concept of robustness to mutations is well established in 
the literature CSlIISl. It is intimately tied to neutral correla¬ 
tions, in fact robustness helps quantify the amount of neutral 
correlation present. Before we study the more novel topic of 
non-neutral correlations, it is therefore interesting to compare 
various measures of robustness between the biological GP maps 
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FIG. 2. Greater mutational robustness indicates the presence of neutral correlations. A) The phenotype robustness pp is plotted as a function 
of frequency fp for all phenotypes in the RNA secondary structure models: RNA12, RNA 15, RNA20, the Polyomino models for protein quaternary 
structure: S 2,8 S 3,8 and the HP protein folding model HP24. Each model has an associated random model with the same frequencies, but we only 
show one example, with K = 4 and L = 12 and a set of phenotypes chosen with a broad range of frequencies to best illustrate the relationship 
(red points). All random models closely follow the expected theoretical curve pp = fp (grey line). The biophysical models exhibit a much larger 
robustness than the random models, which indicates the presence of positive neutral correlations. The red dotted line is 5 (Eq. lU) for K — 4, 
L — 12. If (p > 5) then large neutral networks are expected, which is much more likely for the biophysical models than for the random model. 

B) The average n-robustness defined in Eq. 3 for each of the three biological GP maps, along with the expected values ^ 

for the associated random null models (fiat coloureoTiorizontal lines) is plotted against n. Across all three GP maps, we see a typical decay in 
robustness towards the random null model expectation with increasing mutational distance. Prom this decay a neutral correlation length can be 
defined which is shorter for the HP model than for the other two models. Error bars for HP24 are the standard error on the mean of the average 
n-robustness. 


Properties shared by random and biological GP maps 

Symbol 

Alphabet size: 

K 

Genotype length’. 

L 

Number of 1-mutation neighbours of a genotype: 

{K-1)L 

Number of genotypes: 

No = 

Number of phenotypes: 

Np 

Redundancy, the size of neutral set or the number of 
genotypes that map to phenotype p 

Fp 

Phenotype frequency, the fraction of genotypes that map 
to phenotype p. 

II 

Properties that differ between random and biological GP 
maps 

Symbol 

Neutral set. all genotypes that map to phenotype p 

Qp 

Neutral component A subset of Qp that is fully connected 
by point mutations. Also called a neutral network. 

NN 

The number of 1-mutation neighbours of genotype g 
mapping to phenotype p 

'^P,9 

Phenotype robustness’, mean robustness of all genotypes 
mapping to a phenotype p 

Pp 

Phenotype mutation probability’. Probability that a point 
mutation from a genotype mapping to phenotype p will 
generate a genotype mapping to phenotype q 

4>qp 


TABLE I. GP map nomenclature. 


and the random uncorrelated GP map. 

The 1-robustness of a single genotype g that maps to pheno¬ 
type p is straightforwardly defined as the number of genotypes 
rip^g that map to p that are accessible within one point muta¬ 


tion of g. The phenotypic robustness pp of a phenotype p is 
defined the average of the 1-robustness over the entire neutral 
set Qp 1^ . This can be expressed algebraically as 




Tlnj 




( 1 ) 


In a random GP map, phenotypes are arranged randomly over 
genotypes so the probability that a genotype leads to phenotype 
p is given by its frequency fp, independently of the identity of 
its neighbours. The phenotypic robustness therefore is simply 


Pp — fp 


and the mean number of neutral neighbours is 
(n,,p) = {K- 1)LU 

which is the expectation value for a binomial distribution with 
{K — 1)1/ trials and probability of a given neighbour being fp. 
It is independent of the identity of the genotype g. 

We define neutral correlations as the difference in how geno¬ 
types mapping to the same phenotype are distributed in a bio¬ 
logically relevant GP map as compared to the associated ran¬ 
dom GP map null model. One way of characterising these neu¬ 
tral correlations is by comparing the phenotype robustness pp to 
the random expectation pp = fp. The violation of this equality 
is a sufficient (though not necessary) condition for the existence 
of neutral correlations. Moreover, we define a phenotype p to 
have positive neutral correlations if pp > fp is satisfied. This is 
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intuitive - when robustness is greater than fp then phenotypes 
are closer to each other in the genotype network than would be 
expected by random chance. 

Using the above definitions around neutral correlations, we 
explicitly consider the robustness in the various GP maps. In 
Fig.[^, we compare the phenotypic robustness across our three 
biological GP maps to the robustness of the associated random 
GP map. The figure confirms both the analytical result derived 
above for the random model that pp = fp (we only show one 
schematic random map in the figure, but the others have the 
same behaviour). In sharp contrast, for the biological GP maps 
we find that, very roughly, pp oc log /p, so that the robustness is 
much larger than would be expected for the null model, in fact 
by several orders of magnitude for smaller fp. Since pp ^ fp, 
this indicates the presence of extremely strong neutral corre¬ 
lations in these biological GP maps. Of course the fact that 
Pp > fp is not a new finding, but it is instructive to show this 
trend displayed explicitly for entire mappings in the three kinds 
of biological systems. 


B. Generalised robustness and neutral correlations 


We next extend phenotype robustness to n-mutations. Gen¬ 
eralised robustness or n-robustness p ^^, measures phenotypic 
robustness for a greater number of mutations. It is defined as 
the robustness of a genotype with phenotype pton independent 
mutations to its genotype, rather than just the single mutation 
discussed above. This can be expressed algebraically as 


^(n) _ J_ ^ (n)_ ]_ _ 


( 2 ) 


(n) 

where rip J is the number of n-mutant neighbours of g with phe¬ 
notype p and the normalisation on the right-hand of the sum 
is the total number of n-mutants. In the same way as for the 
phenotype robustness, the n-robustness is averaged across the 
neutral set Qp of all genotypes that map to phenotype p. 

A further quantity we define is the average n-robustness 
which is the average of the n-robustness over all phe¬ 
notypes in a given GP map: 


n-robustness: 



since the phenotype frequencies in a GP map sum to unity. The 
inequality ^ ^/Np can be used to define whether a bio¬ 

logical GP map possesses neutral correlations as a whole. 

We consider the average n-robustness against the radius n for 
the three GP maps 5'2,8, RNA12 and HP24. A sample of 100 
genotypes for each phenotype in the respective systems is taken 
(apart from HP24 where a sample of 100 randomly chosen phe¬ 
notypes is made due to the large number of phenotypes) and 
the n-robustness is measured and averaged over phenotypes. In 
Fig.[^, we plot the average n-robustness at each radius along 
with the fiat expectation lines from Eq. |^for the null models. 
In all three cases we observe a decay from greater than the null 
values for small radii to slightly less than the null expectation at 
larger radii. The reason for this drop below the random expec¬ 
tation can be understood intuitively: given that positive neutral 
correlations are present, the over-representation for small radii 
must be balanced at larger radii by under-representation in order 
for the number of genotypes to balance. 

We also define a neutral correlation length ri" which mea¬ 
sures the mutational hamming distance over which neutral cor¬ 
relations extend. We define n* for a phenotype to be equal to 
the smallest value of n where p^^^ < fp and for the GP map 
when < l/Np. We find that n* = 7 for the RNA12 

model, n* = 6 for the Polyomino 82,8 model and n* = 5 for 
the HP24 model. The neutral correlation length is smaller for 
the HP model than the other two systems. As discussed in the 
methods section, and illustrated in Fig|^, the HP model typ¬ 
ically has phenotypes with smaller frequencies/robustness than 
the other two systems suggesting neutral networks that do not 
expand to the same diameter which would reduce the expected 
neutral correlation length. All three models are of fairly small 
genome length L, so one should be careful of reading too much 
into the numerical values of these correlation lengths. However, 
it may very well be that this ordering of models will persist for 
larger L. 
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Np 
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C. The presence of positive neutral correlations/higher 
(^) phenotype robustness results in larger and fewer neutral 

components 


where V is the set of all Np phenotypes in the GP map. In 
contrast to the two previous definitions that measure robustness 
for a single phenotype, it is a general property of the whole GP 
map. One could imagine generalising this further to a subset 
of the phenotypes, for example those whose frequencies fp are 
greater than the average Np /Nq- 

To establish the n-robustness and average n-robustness in the 
random GP map, the same logic can be applied as in the pre¬ 
vious section. Since the probability of finding a phenotype is 
uniformly distributed over the genotype space, the n-robustness 
is given by 

An) _ r 

Pp — Jp 

with the n-robustness the same for all n, leading to an average 


Having illustrated the concept of positive neutral correlations 
- measured by (generalised) robustness greater than that of the 
random null model - we next show how other properties of neu¬ 
tral networks are affected by their presence. 

The neutral set Qp is the set of all genotypes mapping to phe¬ 
notype p. A component is the subset of the neutral set Qp that 
is connected by single point mutations. We use this term be¬ 
cause it is commonly used in graph theory to denote a set that 
is connected. Although the literature can be somewhat ambigu¬ 
ous, with the term neutral networks sometimes referring to the 
neutral set, and sometimes to a neutral component, we take a 
neutral network to be synonymous to a neutral component in 
this paper because if we have only point mutations then a pop¬ 
ulation can only explore a neutral component and may not be 
able access the whole neutral set. 
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FIG. 3. Biological GP maps have much larger and fewer neutral components than their random counterparts due to neutral correlations. 

A) The logarithm of the largest neutral component for a given phenotype is plotted as a function of frequency for random null models (with A = 4, 
L = 12) and three biological GP maps, RNA12, 82,8 and HP24. The vertical dotted line denotes the giant component threshold S ^ 1/36, defined 
in Eq. ([^, for the schematic random model with A = 4, L = 12. The vertical dashed line denotes the single component threshold A 0.37, 
defined in Eq. for the schematic random model. The biological GP maps show much larger connected components below these thresholds, 
due to the presence of positive neutral correlations. B) The logarithm of the total number of neutral components against frequency is plotted for 
the same models. The theoretical thresholds S and A work well for random model but again the number of components in the biophysical models 
differ greatly from the random model expectation due to the presence of correlations. In both plots, error bars represent a single standard deviation 
from the 100 independent realisations of the random null model used to derive the neutral component statistics. 


There are several reasons why a neutral set may not be fully 
connected by neutral point mutations. If the genotypes are too 
diffusely spread out over the full genotype space, then they 
may be disconnected. But in some cases basic biophysical con¬ 
straints, such as the neutral reciprocal sign epistasis described 
in the Methods section, also lead to fragmentation. 

We begin by comparing the size of neutral components in the 
random null model to those found in our biological GP maps. 
In the random model, there are two important threshold values: 
firstly, the giant component onset, when a phenotype’s compo¬ 
nents change from being largely isolated to forming larger con¬ 
nected clusters, and secondly, the single component onset where 
virtually all genotypes are taken up by a single giant connected 
component. 

As each genotype has many neighbours, a simple mean-field¬ 
like approach from percolation theory for random graphs 
should be fairly accurate. This suggests that the giant compo¬ 
nent onset begins when the average number of neighbours of a 
given genotype with the same phenotype is approximately unity, 
which was also the criterion used by John Maynard Smith (T] . 
For the null model, where phenotypes are assigned to genotypes 
completely randomly, this reduces to an explicit threshold fre¬ 
quency 


{K-1)L 


(5) 


such that we expect the giant components for phenotypes with 
fp > S. It can be shown analytically in the limit L ^ 00 
that there is another transition at 


A = 1 - 


1 

A^ 


( 6 ) 


where, for fp > A, all the components coalesce into one single 
giant component, so that the neutral set should be (nearly) fully 
connected. While the giant component threshold 6 scales as 
1/L, so that it decreases for larger maps, the single component 
threshold A from Eq.j^is independent of genome length L, and 
only varies with alphabet size. For example, A = 0.5 for A = 2 
and A ^ 0.37 for A = 4. These are large frequencies that are 
unlikely to be reached for more than a single phenotype in any 
realistic GP maps. 

In Fig. we plot how the largest component size (left) and 
number of components (right) varies with frequency in both a 
null model (A = 4, L = 12) and three GP maps 82 , 8 , RNA12 
and HP24. We first focus on the simple schematic null model. 
Data is calculated by averaging over 100 independent realisa¬ 
tions of the random mapping of genotypes to phenotypes in a 
way that preserves the frequencies. The largest component size, 
and the number of components formed by the phenotype, are 
then measured. These values are shown in Fig. for an array 
of frequencies in the schematic null GP map. Below the giant 
component onset 6 « 1/36, most genotypes are completely iso¬ 
lated - the total number of neutral components scales with fp. 
Around the giant component threshold 6 , this scaling changes 
markedly, and instead the size of the largest neutral components 
scales linearly with fp and takes up the majority of the geno¬ 
types in the neutral set. The number of components continues 
to decline until fp exceeds the single component connectivity 
threshold A ^ 0.37, at which point there is just one component 
and the neutral set is completely connected. 

We next consider the biological GP maps relative to the be¬ 
haviour exhibited by the null model. Firstly, all three GP maps 
have much larger maximum neutral set sizes than the random 
model. This is not surprising, as Fig.[^ shows that, due to pos- 
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itive neutral correlations, Pp > S for most phenotypes in each 
system (p = 6 for K = 4, L = 12 is shown as a dotted red line 
in the plot). Once the probability of having a neutral neighbour 
is above the S threshold, we expect large networks. For HP24 
and RNA12, the largest neutral component size clearly grows 
linearly with frequency, and so scales linearly with the size of 
the neutral set. For the Polyomino space £' 2,8 this scaling is 
less evident, but the components are still much larger than their 
random counterparts would be. 

Secondly, for all three models, the number of components 
does not vary much with fp, in contrast to the random model 
where this number scales, as expected, with the neutral set size 
if fp ^ Since these components typically have robustness 
above 6 or even A, the reason there are still multiple components 
must be due to biophysical constraints which are not present in 
the random model, such as the neutral non-reciprocal sign epis- 
tasis discussed earlier for RNA. These effects are to first order 
independent of fp which explains why the number of compo¬ 
nents does not correlate with fp. In each of these three models 
the largest “phenotype” of all is the deleterious non-folding or 
non-assembling one. Its frequency exceeds the threshold A and 
its neutral set is fully connected. 

Differences between the three biophysical GP maps observed 
in Fig. can be fairly easily explained by some of the differ¬ 
ences highlighted in the methods section. For example, the 
number of phenotypes per genotype is largest in the HP model, 
and smallest in the Polyomino model, which explains why they 
group at different frequencies. While the number of compo¬ 
nents in Fig. B) is generally much lower in the biological 
models than it is in the random model, the HP model has signifi¬ 
cantly fewer components than the RNA and Polyomino models 
do. This may be due to the fact that the HP model does not 
exhibit effects such as neutral reciprocal sign epistasis, which 
fragments neutral sets in the other two systems. 

We conclude that due to positive neutral correlations and the 
concomitant higher robustness, the biophysical models consid¬ 
ered here have large neutral networks even for frequencies fp 
that are several orders of magnitude lower than the random 
model large component threshold S. The abject failure of the 
random model to predict the robustness and the neutral network 
size highlights the importance of neutral correlations in these 
systems. 

D. Non-neutral phenotype mutation probability 

We next consider non-neutral mutations. The first question is: 
Are two different phenotypes, on average, more or less likely to 
be connected to each other than one would expect by chance? 
To address this question, we employ a generalisation of robust¬ 
ness, namely the phenotype mutation probability <fqp of q with 
respect to p, defined as the fraction of 1-point mutations of 
genotypes in the neutral set for phenotype p that map to phe¬ 
notype q. This can be written as: 

^ FJK - l)L ^ 

Thus <fqp averages a local property, riq^g - the number of geno¬ 
types that map to phenotype q found the 1-mutation neighbour¬ 
hood of a genotype that maps to phenotype p - over the entire 
neutral set Qp. Note that this phenotype mutation probability 
is not symmetric (0^^ (fpq) and that, if p = q, it reduces 


to the phenotype robustness (fpp = Pp. It has recently been 
shown ca that (l)qp is a key quantity for incorporating the struc¬ 
ture of a GP map into population genetic calculations. 

In the null model we expect <fqp = fq to be an excellent ap¬ 
proximation ca, with the caveat that it must be possible for 
enough genotypes to be sampled. What do we mean by enough 
genotypes? Given a phenotype p with redundancy Fp, there are 
at most Fp{K — 1)L unique neighbours available. This num¬ 
ber provides an upper bound - in reality, several neighbours of 
one genotype will also be neighbours of another genotype with 
the same phenotype, resulting in a reduction in the number of 
unique neighbours. Nevertheless, this allows us to define a min¬ 
imum threshold 


' Fp{K -1)L 

If/, < 7, then the expected number of genotypes with phe¬ 
notype q found around phenotype p is less than one, and the 
probability that <fqp = 0 due to statistical fiuctuations becomes 
appreciable. Further detail on how <fqp and Fp relate when the 
threshold is not satisfied, which is mainly relevant for smaller 
GP maps and for lower Fp, is provided in Appendix [A[ Here we 
focus on phenotypes with larger Fp, in the larger GP maps of 
the previous section. 

In Fig. 1^ we plot the relationship between the phenotype 
mutation probability (fqp and global frequency fq around the 
RNA20 phenotype with the second largest neutral set, the as¬ 
sembling phenotype for with the largest neutral set, and the 
HP5x5 folding phenotype with the largest neutral set. For phe¬ 
notypes in and HP5x5, with such large numbers of geno¬ 
types, all phenotypes have fq values that are significantly above 
fq = 1 (vertical dotted lines), which is the approximate thresh¬ 
old at which at least one genotype of phenotype q would be 
expected to be found. A small fraction of phenotypes lie close 
to the fq=y threshold for RNA20, but by far the majority may 
be expected to be effectively sampled. For RNA20 and Ss^s, 
we observe a very strong and highly significant positive corre¬ 
lation with the random null model expectation (fqp = fq. In 
HP5x5, there is also a strong positive correlation, though less 
strong than in the RNA and Polyomino cases, with a greater 
number of phenotypes falling below the one-to-one expectation. 
We did not plot the non-compact model HP24 because most of 
its frequencies are below the threshold 7 (see supporting infor¬ 
mation). 

To summarise, in contrast to the robustness pp = f>pp where 
neutral mutations lead to strong deviations from the null model, 
the non-neutral phenotype mutation probabilities follow the 
random model expectation that (j)qp « fq remarkably well. 
There are still important deviations, especially for those pheno¬ 
types that can not be reached due to biophysical constraints so 
that (fqp = 0 1^ . Moreover, it may be an interesting exercise to 
look more closely at phenotypes for which is significantly 
greater or less than fq as such deviations could signal similari¬ 
ties or differences between phenotypes. For example, two RNA 
phenotypes with similar hairpin topology, but perhaps a differ¬ 
ence of one bond in a stem may have a larger probability of 
interconverting than topologically more dissimilar RNA pheno¬ 
types. The difference between (j)qp and fq could then be used 
to quantify the difference between phenotypes p and q. These 
more subtle types of correlation are beyond the scope of this 
paper. At any rate, compared to the result in the previous sec¬ 
tions showing the strength of neutral correlations, the dominant 
agreement with the random model is apparent. However, given 
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FIG. 4. Phenotype mutation probabilities scale with global frequency. We present results for the three GP maps: A) RNA20, B) S 3,8 and C) 
HP5x5. We plot the relationship between (l)qp (circles) and fq for the largest non-deleterious phenotype p in 5'3,8 and HP5x5, and for the second 
largest in RNA20 (not the largest due to computational expense). We see in each case a strong positive correlation (p-value <C 0.05 in all cases), 
very similar to the expectation for the null model (not shown here, but for which the correlation is exact to within statistical fluctuations, see 
ref. |[T4l and Appendix [B]). Spearman rank correlation coefficients are shown in the top-left of each plot. Differences from = fq are relatively 
small compared to the overall range of variation, except for sets of phenotypes that are not connected at all, which typically arise due to biophysical 
constraints. These are shown as downward triangles along the lower horizontal dotted line which represents (fqp — 0. For each plot, the upward 
triangle indicates (fpp — pp, the phenotype robustness, which is always over-represented {pp ^ fp) due to neutral correlations. 


that <f>qp is averaged over a neutral set, it may be that there are lo¬ 
cal non-neutral correlations that are obscured by the averaging. 
With this in mind, we next investigate such local correlations. 


E. Non-neutral local over-representation correlations 

We first describe non-neutral local over-representation cor¬ 
relations which mean that, given phenotype q is found in the 
1-mutation neighbourhood of a genotype g (which maps to phe¬ 
notype p q), then phenotype q will appear a greater number 
of times in total than predicted by fq or <f>qp in this 1-mutation 
neighbourhood, as pointed out in ref. 12^ . These correlations 
are illustrated in Fig.[^. 

To measure 1-mutation neighbourhoods, we sample ran¬ 
domly chosen genotypes g from the neutral set Qp, with a geno¬ 
type of phenotype q in its neighbourhood. We then measure the 
phenotype of all other neighbours of g. From this sample, we 
obtain the probability P{q,p^ m) of q occurring m times in the 
1-mutation neighbourhood of a genotype mapping to phenotype 
p, given that q occurs at least once. 

Two control null expectations may also be derived for 
P{q,p,m). In the random model where phenotypes are ran¬ 
domly assigned, given q is in the 1-mutation neighbourhood of 
a genotype g (at a specific genotype g'), the probability may 
be calculated as a binomial probability based upon the overall 
frequency of q, leading to 

A second null expectation calculates the binomial probabil¬ 
ity by replacing fq in Eq. above by replacing the phenotype 


A) Local over-representation of phenotypes 


null model 

• • 



biological GP map 



B) Mutational neighbourhood correlations 

neutral neighbours 




neutral but not 
neighbouring 


FIG. 5. Illustration of further non-neutral correlations. A) On 

the right, the orange phenotype is over-represented relative to the null 
model: The red genotype in the centre has more orange neighbours 
than would be expected by the global frequency of orange. B) The 
phenotypes that appear in the mutational neighbourhood of two neu¬ 
tral neighbours are expected to be more similar (right) than two non¬ 
neighbouring genotypes of the same phenotype (left). 


mutation probability (fqp for the GP map instead: 

P 2 {q,P,m) = 

(9) 

In contrast to Pi{q^p,m), this form accounts for any overall 
phenotypic heterogeneity known to be present in the GP map. 
We compare the actual local prevalence against these two null 
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expectations in Fig. For RNA20, Ss^s and HP5x5 we chose 
the same three phenotypes for phenotype q as we did in the 
previous section, while for phenotype p we choose one-by-one 
the next n = 10 largest (non-deleterious) phenotypes available 
in the GP map. By sampling 10,000 neighbourhoods for each 
of the n = 10 phenotypes for p, we calculate an average for 
P{q, p, m) across the phenotypes (P{q^ p, m)) and compare this 
in Fig. ^ to the averages for the null expectations Pi{q,p^m) 
and p 2 (q^ p, m ). For each biological GP map, q is more likely to 
be over-represented, that is to appear multiple times if it appears 
at least once when compared to the null expectations, leading to 
a skewed distribution compared to the control case. The most 
striking result is seen in RNA20, where there is a substantial tail 
to the distribution. We use average measures here to provide 
the general profile, smoothing out particular features that may 
occur between individual pairs of phenotypes q and p, but the 
local over-representation is seen for any of the phenotype pairs 
considered. 

One consequence of these local over-representation correla¬ 
tions is that the probability that a genotype with phenotype p has 
phenotype q at least once in its 1-mutant neighbourhood is less 
than expected from (pqp. This is because those genotypes that 
have phenotype p in their 1-mutant neighbourhood typically do 
so a greater number of times than expected. This must therefore 
be compensated with fewer genotypes around which phenotype 
p actually appears at all, which we confirm numerically. In the 
RNA20 GP map, with the most frequent of the set of phenotypes 
used for p and the next most frequent used as g, the probability 
of finding q at least once is 0.12 versus a null expectation of 
1 — (1 — = 0.20. Thus these correlations lead to 

heterogeneity in the connections between phenotypes. 

How these correlations affect evolutionary dynamics will de¬ 
pend on the regime being explored O. If the population is 
neutrally exploring genotypes that map to phenotype p, then 
in the monomorphic regime of evolutionary dynamics, where 
NLji <C 1, this heterogeneity will lead to a significant drop 
in the rate at which q is first discovered by neutral exploration. 
In the polymorphic regime where NLji ^ 1, and different in¬ 
dividuals in the population have different genotypes, the rate 
at which novel variation with phenotype q occurs may not be 
that different from the expectation given by (pqp, at least if the 
population is spread across a large enough number of different 
genotypes to average over local heterogeneity lIT^ . 


F. Non-neutral local mutational neighbourhood correlations 

We next examine non-neutral local mutational neigh¬ 
bourhood correlations which are illustrated schematically in 
Fig.|^. They show that the 1-mutation neighbourhoods of two 
genotypes connected by a neutral point mutation are more likely 
to have similar phenotypic compositions than would be ex¬ 
pected by two randomly chosen neutral non-neighbouring geno¬ 
types of the same phenotype. This type of correlation has al¬ 
ready been demonstrated to exist for RNA |[39ll . To measure the 
similarity of neighbouring genotypes’ mutational neighbour¬ 
hoods, we consider the local quantity {K — 1)1/, 

which becomes (pqp when averaged over the whole neutral set 
Qp. We compare the for neighbouring genotypes with 

non-neighbouring genotypes in both the null model and biolog¬ 
ical GP maps. The similarity or difference could be measured 
in several different ways. The statistical measure we employ 


here is the Bhattacharyya coefficient PlOl . which for two dis¬ 
crete probability distributions Xi and yi may be expressed as 

BC{xi,yi) = (10) 

i 

varying between 0 and 1 for maximally dissimilar and identical 
discrete probability distributions respectively. 

To quantify whether neutral neighbours g and h have 
more similar phenotype distributions in comparison to non¬ 
neighbouring neutral genotype pairs g and g2, we com¬ 
pared the similarity ratio of the Bhattacharyya coefficients, 
BC{g,h)/BC{g, g2), using the to define the distribu¬ 

tions. A ratio greater than unity indicates that the phenotype 
distributions around neutral neighbours are more similar than 
the randomly selected neutral pair, and vice versa. We remove 
the AT — 2 mutual neighbours of g and h from the distributions 
as these will automatically contribute to similarity between the 
neighbourhoods in a trivial manner which we wish to exclude. 

In Fig. [7] we plot histograms of the similarity ratio for 10,000 
samples of g,h and g2 in RNA20, and HP5x5, where the 
phenotype sampled has the second largest frequency in RNA20, 
and the largest frequency in Ss^s and HP5x5 (excluding the del 
phenotype). For 10,000 samples the means are 1.357±0.003 for 
RNA20, 1.063 ± 0.001 for 5'3,8 and 1.025 ± 0.001 for HP5x5, 
where the error is the standard error on the mean. For RNA20 
and 5'3,8, a clear skew in the overall distribution may be visually 
observed, demonstrating that neutral neighbours, on average, 
have more similar mutational neighbourhoods. HP5x5 also has 
the mean of its distribution at a value slightly larger than unity 
but it is much more marginal in this case, and the skew is harder 
to detect. We note that in general, the non-neutral correlations 
are weakest for the HP5x5 model. Finally, just as is the case 
for the non-neutral local over-representation correlations of the 
previous section, these local mutational neighbourhood correla¬ 
tions also reduce the rate at which novel phenotypes would be 
discovered by neutral exploration since a neutral neighbour is 
more likely to have some of the same phenotypes in its muta¬ 
tional neighbourhood, and so fewer alternatives. 


G. Non-neutral deleterious phenotype correlations 

The final, and perhaps most important, type of non-neutral 
correlation we consider is the accessibility of the deleterious 
phenotype from folding or self-assembling phenotypes, which 
we call non-neutral deleterious phenotype correlations. This 
type of non-neutral correlation is closest to the type of correla¬ 
tion suggested by Maynard Smith m 

In Fig. we plot histograms of the ratio 0dei,p//dei for all 
phenotypes p in and HP5x5, and the top 20 most frequent 
(largest fp) in RNA20 (limited due to computational expense of 
this larger system). In all cases, we see that the deleterious phe¬ 
notype is significantly less frequent around the non-deleterious 
phenotypes. This behaviour contrasts to non-deleterious phe¬ 
notypes, for which <pqp ^ fq. As sl corollary of this effect, we 
also find Pdei//dei equal to 1.10,1.16 and 1.19 for RNA12, Sg^g 
and HP5x5 respectively, illustrating positive correlations, that 
is, a corresponding local over-representation of the deleterious 
phenotype in its own mutational neighbourhoods. Moreover, 
for 1/ = 20 we find that Pdei//dei = 2.34 suggesting that these 
positive neutral correlations may become stronger for larger L. 
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FIG. 6. Non-neutral local over-representation correlations result in phenotypes being more likely to be found multiple times around 
genotypes. We present results for the three GP maps: A) RNA20, B) 83,8 and C) HP5x5. We pick the same frequent phenotypes q in each of our 
biological GP maps as used in Fig.|^ and consider the prevalence of q around genotype g with phenotype p, given that q occurs at least once in 
the 1-mutation neighbourhood of g. The average of P{q,p, m) across the n = 10 most frequent phenotypes p in the neighbourhood of q (with 
p ^ q and p / del), is compared to the respective averages for random null expectations Pi{q,p, m) and p 2 {q,p, m) defined in the text. The 
mean of each distribution is plotted as a dotted line in each case. Contiguous sections with a probability greater than 10“^ are joined with lines 
in order to guide the eye. The mean value of m for each of the biological GP maps and the two random controls are shown as respective dotted 
lines with the same colours. Compared to the two null expectations of occurrence, q is over-represented locally as demonstrated by the shift of the 
means to the right. 



Similarity ratio, 



FIG. 7. Non-neutral local mutational neighbourhood correlations result in mutational neighbourhoods of neutral neighbours being more 
similar than randomly selected neutral pairs. We present results for the three GP maps: A) RNA20, B) 83,8 and C) HP5x5. Using the ratio of 
Bhattacharyya coefficients defined in Eq. we show that neutral neighbours (p and h) have a closer phenotype probability distribution than a 
randomly chosen neutral pair (g and P 2 ). This is seen through the ratio being skewed with a mean (coloured vertical dashed lines) larger than unity 
(black vertical dashed lines). The standard error on this mean is negligible compared to the distance of the mean from one. 


We find that the del phenotype forms only a single compo- nent even in the random null model, 
nent in RNAI2, 82 , 8 , HP24 and HP5x5. This result is unsur¬ 
prising because the large size of the del phenotype in each GP 
map (85%, 54%, 98% and 82% respectively) means that the fre¬ 
quencies are all well above the single component threshold A of 
Eq. which would lead to the expectation of a single compo- 

































12 



Deleterious phenotype ratio, 


0del,p 

/del 


FIG. 8. Non-neutral deleterious phenotype correlations: The deleterious phenotype is under-represented in the neighbourhood of folding 
or self-assembling phenotypes. We present results for the three GP maps: A) RNA20, B) Ss.s and C) HP5x5. Histograms of the ratio of the 
phenotype mutation probability (0dei,p) divided by the null model expectation of the global frequency (/del) for the deleterious phenotype (non¬ 
folding for RNA/HP, non-assembling for Polyominoes). The distribution is clearly skewed to values < 1, as highlighted by the dashed vertical 
coloured lines representing the mean in each case. 


V. DISCUSSION 

In this paper, we have explored the role of genetic corre¬ 
lations, which we defined and quantified as the difference in 
how genotypes are mutationally connected for biologically rel¬ 
evant GP maps, compared to a random null model with the same 
global properties (alphabet size, genome length, and number of 
genotypes per phenotype). Genetic correlations provide a sim¬ 
ple conceptual framework within which a number of topological 
properties of GP maps can naturally be captured. 

Neutral correlations: We first explored the phenotype ro¬ 
bustness for all three GP maps, showing that pp > fp for all 
phenotypes, a result which is not unexpected in the literature, 
but to our knowledge has not been compared for a set of whole 
GP maps before. Since pp = fp for the random model, the ex¬ 
tent to which Pp differs from fp can be viewed as a measure of 
the extent of the neutral correlations. 

We also introduced the concept of n-robustness, which mea¬ 
sures robustness over n mutations. From this we derived an¬ 
other criterion that measures the presence of neutral correlations 
by averaging this measure over all phenotypes and comparing 
to the null expectation: If > 1/A/p then there are posi¬ 

tive neutral correlations. We find that the enhanced probability 
of encountering a genotype mapping to the same phenotype can 
extend to multiple mutations n away from genotypes. The ex¬ 
tent of the correlations in sequence space can be quantified by 
the number of mutations n* at which the criterion is violated, a 
measure we call the correlation length of the neutral mutations. 
We find that that n* is largest for the RNA12 model, and small¬ 
est for the HP24 model. How n* or even the relative ordering of 
the correlation lengths between the different systems will scale 
with increasing genome length L remains an open question. 

It should also be emphasised that the full complexity of neu¬ 
tral correlations for a phenotype are only partly be captured with 
the measures we introduced here, which average over the entire 


neutral set. As can be seen in Fig 3A of ref. ifT^ . a single neu¬ 
tral set can have significant local heterogeneity in its internal 
connections. Since neutral sets are frequently so vast that they 
cannot be fully explored by populations on evolutionary time- 
scales, such local heterogeneities may also have implications 
for evolutionary dynamics. Thus local measures of heterogene¬ 
ity and robustness, which may be infiuenced, for example, by 
the identities of multiple neighbours, or the position of a muta¬ 
tion along a genome, may also be important to develop in future 
work. 

We found, for the three biological GP maps that we studied, 
that the dominant relationship of robustness with frequency is 
Pp ^ log fp, a scaling that has already been pointed out earlier 
for RNA |[T3][4T1. In an interesting paper that applies concepts 
from network theory |[38]| to neutral sets, Aguirre et. al. na ra¬ 
tionalise this scaling for RNA by separating out the mutational 
behaviour of bound and unbound bases. It would be interesting 
to see if a more general argument could be developed to ex¬ 
plain the logarithmic scaling across all the systems we studied. 
Moreover, these results also pose fascinating questions relating 
to why or how the constraints in a GP map lead to the kinds of 
neutral correlations they do. 

Some clues to the underlying causes of neutral correlations 
and robustness can be gleaned from Maynard Smith’s origi¬ 
nal paper, where he illustrated the concept of a neutral net¬ 
work with the parlour game of transitioning between connect¬ 
ing two words in the English language by changing one letter 
at a time, with each change also generating a valid word. He 
used the example of chan^g “WORD” to “GENE” in four 
steps as illustrated in Eig.^ There are 26^ = 456,976 dif¬ 
ferent possible 4-letter words, but, according to the Merriam- 
Webster Official Scrabble Players Dictionary 2014, only 4,175, 
or just over 0.9%, are valid English words. If we consider the 
set of valid words to be a phenotype, then it has frequency of 
only fp = 0.009, just under the giant component threshold 
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WARD 


GENT 


WOLD 


GENU 



GENS 


GORE- 


GONE- GENE 


GANE 


BENE 


DENE 


SENE 


NENE 


FIG. 9. Enhanced robustness in Maynard Smith’s 4-letter word game. The single mutation path WORD ^ WORE ^ GORE ^ GONE ^ 
GENE is shown in red. All valid words within a one letter mutation of “WORD” and “GENE” are also depicted. According to the Merriam-Webster 
Official Scrabble Players Dictionary 2014, only 4,175 of the 456,976 possible 4-letter words are valid English words (at least for Scrabble). Since 
each word has 100 neighbours, for a random model, the expected number of valid words within a one letter mutation is < 1. Nevertheless, due 
to positive neutral correlations, the probability that a valid word has another valid word as a 1-mutant neighbour is more than ten times greater, 
as illustrated above for “WORD” and “GENE”. Such correlations make the game much easier. As pointed out by Maynard Smith, in a biological 
system, such correlations (in his case between “meaningful sequences”) can facilitate evolutionary dynamics m 


5 = l/(Ar — 1)I/ = 0.01, and well below the single component 
threshold A ^ 0.12. On average the probability that a 4-letter 
word has a valid neighbour is just below one. However, if we 
measure the phenotype robustness, that is the mean probability 
that a valid 4-letter word has valid words in its (FT — 1)1/ = 100 
neighbours, we find that pp ^ 0.11, or on average each word 
has 11 neighbours, which makes the game much simpler than if 
the random expectation held. 

The reasons this game exhibits such a large enhancement of 
the robustness over fp clearly arise from neutral correlations in 
English words. For example, vowels are more likely to appear 
in specific places in words than would be expected by chance. 
The second letter of 4-letter English words has a 74% chance of 
being a vowel compared to the 5/26 = 19% overall average prob¬ 
ability per locus. So if a word has a vowel placed at the second 
letter, it is much more likely to have neighbouring words using 
the same vowel, as can be seen in Fig. This example illus¬ 
trates how basic properties (in this case of language properties, 
in the case of our models, biophysical properties) can generate 
correlations, which can then result in high robustness. 

A question often debated in the literature is the extent to 
which mutational robustness is selected for. Here we argue that 
a major enhancement of robustness, often by many orders of 
magnitude over the random expectation, is not caused by selec¬ 
tion, but rather emerges from the internal constraints of a GP 
map - the way that genotypes map to phenotypes - which natu¬ 
rally lead to positive neutral correlations. It may still be the case 
that more robust genotypes can be selected for within a neu¬ 
tral set, or that these genotypes are favoured in certain dynamic 
regimes 14^ . It may also be true that in some cases a particu¬ 
lar phenotype is preferred by selection because it is more robust 
than an alternative one. But even if this is so, it is important 
to keep in mind that natural selection is still acting on variation 
that is already naturally quite robust due to correlations caused 
by biophysical constraints. 

The relationship between robustness and evolvability has 
been the subject of much discussion in the literature (121 ESI 
l43l l44l . Here we show, as already anticipated by Maynard 
Smith Cl, that if the phenotype robustness is roughly larger than 

6 = 1/{K — 1)L, so that the expected number of neutral neigh¬ 
bours is greater than one, then the phenotype will exhibit large 
neutral networks. In the random model, large networks will typ¬ 


ically be very rare, but neutral correlations mean that robustness 
above the 6 threshold is common for the biophysical GP maps. 
The effect can be very large. For example, for L = 55 RNA, a 
recent study ll33l suggests that there are about Np ^ 8 x 10^^ 
phenotypes, so that the mean frequency is fp ^ 10“^^. In 
fact all phenotype frequencies are well below the threshold 
^ = 1/(3 X 55) = 0.00606 above which we expect extended 
neutral networks. On the other hand, the mean robustness of all 
phenotypes was estimated to be ^ 0.14 > 6 ^ fp. Neu¬ 
tral correlations increase the probability of a nearest neighbour 
generating the same phenotype by on average about 12 orders of 
magnitude over the mean expectation of the null model, lifting 
robustness well above the threshold 6. Thus the most important 
way that neutral correlations contribute to evolvability is by nat¬ 
urally creating robustness greater than the threshold needed to 
generate percolating networks which provide access to pheno¬ 
typic novelty. In fact it may very well be that without neutral 
correlations and its attendant robustness, evolution as we know 
it would not be possible 

Non-neutral correlations: Non-neutral mutations are im¬ 
portant for the generation of novel variation. For all three GP 
maps, the probability 4>qp that a phenotype q is found by a point 
mutation from genotypes mapping to phenotype p is, to first or¬ 
der, given simply by the global frequency: (fqp ^ fq, which is 
independent of p. 

Since fq can span many orders of magnitude, the rate at 
which variation appears (which scales as ^ 1 /(jyqp ^ \/fq 
if a population is neutrally exploring phenotype p (141) can 
also range over many orders of magnitude in these systems. 
These large differences can lead, both in the monomorphic and 
polymorphic regimes, to effects such as the arrival of the fre¬ 
quent na, where frequent phenotypes (with larger fq) fix in a 
population even when alternate phenotypes that are much more 
fit, but much less frequent, are accessible in principle. 

The reason these fitter phenotypes are not fixed is because 
they are unlikely to be found on evolutionary time-scales. Nat¬ 
ural selection can only work on variation that actually arises. 
In the alternative case where the system is effectively in steady 
state, so that a less frequent phenotype has a realistic probabil¬ 
ity to arise in a population, it can still be the case, especially 
at larger mutation rates, that a phenotype with lower fitness but 
larger frequency (and robustness) will fix, an effect known as 






14 


the survival of the flattest na. 

Finally, we note that fqp can be viewed as a non-neutral gen¬ 
eralisation of the phenotypic robustness, but that fpp = pp 
scales very differently with fp than fqp does when p q. In 
the latter case local correlations more or less cancel out when 
averaged over the whole neutral set, so that fqp ^ fq, while 
in the former case the local correlations do not cancel out at all 
because robustness is fundamentally a local quantity. 

It is quite striking that in all three models, a very large num¬ 
ber of phenotypes are indeed connected to one another. The HP 
model merits further discussion in this regard. In a recent re¬ 
view IH, RNA space was compared to “a bowl of spaghetti”, 
because the neutral spaces were connected to most other phe¬ 
notypes, while proteins were compared to a“plum pudding”, 
where the neutral networks were more likely to be isolated from 
one another. We indeed find that the neutral networks in the 
HP24 model are not well connected, but locate the origin of 
this effect in the large NpfNc ratio for HP24, which means 
that many networks are below the threshold of Eq. 0 for con¬ 
nections. By contrast, the compact HP5x5 model with many 
fewer phenotypes but a similar sized genotype space is well con¬ 
nected, more like “spaghetti” than like a “plum pudding”. What 
happens for real proteins, without the simplifying assumptions 
and small system sizes typically studied in the HP model ll34t . 
remains an open question. 

Another type of heterogeneity in the mapping of genotypes 
to phenotypes can be quantified as local non-neutral correla¬ 
tions, which occur when the local neighbourhood of genotypes 
are different from the global expectation given by fqp or fq. 
We investigated two types of correlation (although one could 
imagine many more): i) non-neutral local over-representation 
correlations which result in phenotypes being more likely to 
be found multiple times around genotypes, and ii) non-neutral 
local mutational neighbourhood correlations, which mean that 
two genotypes connected by a neutral point mutation have mu¬ 
tational neighbourhoods that are more similar than do two ran¬ 
domly selected genotypes in a neutral set. 

These two types of correlation mean that the diversity of phe¬ 
notypes in the direct neighbourhood of a genotype is lower than 
expected from the random model or even from the averaged 
phenotype mutation coefficients fqp. Thus the rate at which 
a neutrally exploring population encounters novel variation will 
be reduced due to these correlations. How this effect influences 
evolvability is complex, because the term is used in many differ¬ 
ent ways in the literature ESIIItHID. One type of evolvability 
simply measures the number of different phenotypes that are 
connected by single mutations to a neutral set ll22l . While non¬ 
neutral correlations may not affect this number very much, they 
will affect the rate at which neutral exploration finds these new 
phenotypes. This lowering of the rate at which novelty appears 
may have a larger impact on other measures of evolvability. 

Each of the three models has a deleterious phenotype which 
either does not fold (for RNA and the HP protein model) or 
does not properly assemble (in the Polyomino model for pro¬ 
tein clusters). The third type of non-neutral correlations we 
considered were iii) non-neutral deleterious phenotype corre¬ 
lations. Eor all three GP maps, the folding or assembling phe¬ 
notypes have fewer mutational connections to the deleterious 
phenotypes than would be expected by the global frequency 
/del. This last result is perhaps the most interesting type of non¬ 
neutral correlation. It was already predicted by John Maynard 
Smith in his classic 1970 paper [H, where he argued that ’’mean¬ 


ingful" proteins were more likely to be neighbours of other 
"meaningful" proteins, and by extension, that the probability 
of finding a deleterious phenotype in the mutational neighbour¬ 
hood of a "meaningful" protein would be less than by random 
chance. Such an effect can enhance evolutionary dynamics, be¬ 
cause non-deleterious phenotypes are more strongly connected 
by mutations than expected by random chance, and so the popu¬ 
lation can more easily access potentially meaningful novel vari¬ 
ation. Of course in practice, whether or not even the folding 
or self-assembling phenotypes are in fact “meaningful” will de¬ 
pend on the environment and other factors, but to first order a 
reduced propensity to mutate to manifestly deleterious pheno¬ 
types should be an evolutionary advantage. 

Evolvability: While the effect of neutral correlations on ro¬ 
bustness is straightforward, how correlations affect evolvability 
is more complex, not just because the concept itself is more 
diffuse, but also because the relationships between correlations 
and evolvability are more varied. Nevertheless, we can sum¬ 
marise how different correlations affect evolvability as follows: 

1. Positive neutral correlations, measured by the presence 
of greater phenotypic robustness than would be expected 
by chance, is critical for the formation of large neutral 
networks. These networks are, in turn, a key facilitating 
factor for the ability of a population to access novel vari¬ 
ation by neutral evolution over the network ESIITIIII]- 

Without positive neutral correlations, and the 
associated phenotype robustness, evolvability would be 
hugely suppressed. 

2. The non-folding or non-assembling deleterious set of 
phenotypes are (positively) neutrally correlated and anti¬ 
correlated with the set of folding or assembling pheno¬ 
types. This correlation increases the potential phenotypic 
variation that is accessible by reducing the likelihood of 
mutations leading to a seriously deleterious phenotype. 

3. Local non-neutral correlations generally mean that the 
amount of novel variation available after mutations is 
smaller than one might expect from a random model. 
These correlations will reduce evolvability. Eor exam¬ 
ple, Huynen et al. 1521 showed that the innovation rate 
for a random walk on the neutral space of the L = 76 
secondary structure of tRNA^^® is about 20, even though 
each genotype has 31/ = 228 neighbours. This significant 
lowering of the innovation rate is due to local non-neutral 
correlations, and may be a more general effect. 

4. Eor all three GP maps we found that fqp ^ fp, which 
means that the exponential variation in the frequency of 
phenotypes is refiected in the probability that a novel phe¬ 
notype is found by mutations. As argued in ref. O, for 
a range of evolutionary scenarios, phenotypes with large 
fp will be exponentially more easy to find than pheno¬ 
types with smaller fp. These large differences in the rate 
at which novel variation arrives will strongly modulate 
the evolvability. 

Other measures of genetic correlations: While we have 
built upon and introduced a number of metrics for genetic corre¬ 
lations in GP maps, the framework of correlations has other av¬ 
enues for future computational work. Eor example, the central 
focus here has been on the way phenotypes are neutrally con¬ 
nected and non-neutrally connected without any broader con- 
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cern for the properties of the phenotypes themselves. Pheno¬ 
types themselves have measurable properties such as symmetry, 
size and modularity and one could take the analysis further by 
considering whether there is a relationship between the muta¬ 
tional connections between two phenotypes and their similarity 
based on such properties. Making use of the null model again, 
where there is no predisposition for which phenotypes are mu- 
tationally close together, such phenotype similarity correlations 
could be studied in the biological GP maps. 

Measuring correlations in experiment: Computational 
studies on theoretical systems ultimately need to be backed up 
with empirical evidence in real biological systems. Robustness 
to mutations in protein tertiary structure has been a well-studied 
area in this regard, with both mutagenesis and phylogenetic ex¬ 
periments being used to illustrate robustness ca. It may be 
that non-neutral local correlations could be verified using mu¬ 
tagenesis experiments. For example, for mutational neighbour¬ 
hood correlations, two neighbouring genotypes both with a cho¬ 
sen structure could have their neighbourhoods examined for the 
range of phenotypes and compared to the neighbourhood of a 
more distant genotype with the same methodology used here 
computationally, potentially replicating the findings Fig. but 
for real molecules. It may also be possible to measure the neu¬ 
tral correlation length n* by doing multiple mutation experi¬ 
ments.. However, because it is hard in practice to extract full 
GP map properties such as fp for a given phenotype, the most 
challenging aspect of such an experiment would be in generat¬ 
ing an appropriate effectively random expectation. This same 
challenge holds for the other kinds of correlations we investi¬ 


gate in this paper. 

A few final caveats are in order. In these models it is nat¬ 
ural to use a restricted definition of a neutral mutation leading 
to exactly the same phenotype, whereas a more complete the¬ 
ory would count all mutations that are not visible to selection as 
effectively neutral. Thus the full picture of how these correla¬ 
tion affect evolutionary dynamics is complex, and depends not 
just on the GP map itself, but more generally on the genotype 
to phenotype to fitness map, for which the environment plays 
a key role. Moreover, population genetic parameters such as 
the population size and the mutation rate must be taken into ac¬ 
count. But notwithstanding these complications, the important 
influence that structure in the GP map, in this case measured 
through the lens of genetic correlations, has on the manner in 
which variation arises (the “arrival of the fittest" US), and so 
on evolutionary dynamics, should be evident, confirming May¬ 
nard Smith’s suggestions from many years ago. It may even be 
that without these correlations, Darwinian evolution, and there¬ 
fore life itself, would not have been possible. 
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Appendix A: Associated analytical results of the sampling 
threshold 


We begin with an analysis of the sampling threshold 7 from 
average phenotypes in a given GP map. We showed in Eq. 0 
that for effective sampling, we require 

1 

^ Fp{K- 1)L 
which may be expressed alternatively as 

hh > kl{K-1)L 


The average phenotype frequency may be written down as 


(/) 


phenotype 




1 


or with respect to the genotype sampling distribution as 


(/) 


genotype 




Substituting in the smaller of the two, the average phenotype 
frequency, and then considering the required threshold for ef¬ 
fectively sampling phenotype q from the average phenotype p, 
we find 


^ K^{K - l)L 

For RNA, where empirical scaling values are known {Np ^ 
1.5 X i- 2 i. 8 ^ mi), we can further write 

/,> 0 . 45 ^^ 

L‘2 

for a phenotype q to be effectively sampled. 

We can change the question of effective sampling to ask the 
conditions on Np for the average phenotype q to be accessed 
from the average phenotype p. In these circumstances, we can 
see that 


Np < ^K^{K-1)L 

for which we can see RNA satisfies to an increasing extent for 
increasing L, as 2^31/ — 1.93^ monotonically increases with 
increasing L. 
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Finally, we can also write down the approximate fraction less 
than the average phenotype frequency that is accessible from 
an average phenotype, through expressing fg = x (/)phenotype’ 
which leads to a threshold fraction 


^ ~ K^{K -1)L 


(Al) 


Again, using RNA as an example system this leads to 


X oc 0.81^ 


1 

lA 


for a given length of RNA, showing that an increasing fraction 
of phenotypes with frequencies below the average may be ef¬ 
fectively sampled from the average phenotype. 


Appendix B: Extended analysis of (pqp across a broader range of 
phenotypes 

In this section, in contrast to the analysis in the main body 
considering the frequency of phenotypes q around phenotype 
p where fp 7 , we consider different p across a range of 
phenotype frequencies fp in the GP map. In the random null 
model, at values of fq > 7 , we expect phenotypes to almost 
exactly follow <fqp = fq. When fq < 7 , there are two likely 
possibilities for a given phenotype: either the phenotype is not 
found at all (cjyqp = 0) or it is found a single time (0^^ = 7 ). The 
latter is an over-representation of the local prevalence of q for 
the GP map, while the former is clearly no local representation 
at all. 

In Fig. 10, we display three pairs of plots for the RNA 12, 
£' 2,8 and HP24 GP maps and a randomised null model coun¬ 
terpart. The null models are displayed on the left, with the 
actual GP maps on the right. In each plot, we show the val¬ 
ues of (fqp against fq for three different phenotypes p (coloured 
by data point as red, blue and green, with red the largest fre¬ 
quency phenotype and blue the smallest). Upward triangu¬ 
lar data points represent values for <fpp, downward triangular 
data points <fqp = 0 (shown at le -8 for visualisation purposes 
only) and the circular data points are all other phenotypes. Ver¬ 
tical and horizontal dotted coloured lines represent fq = J 
and <fqp = 7 respectively. The diagonal dashed black line is 
(fqp = fq, the null expectation for phenotypes with fq > 7- 


We begin by discussing the behaviour of the null model. The 
£' 2,8 and HP24 null models provide the extreme cases. For 82 , 8 , 
all phenotypes are highly frequent and have fq Conse¬ 

quently, we see that each of the three phenotypes follows the 
expected trend of f>qp = /^ to a very high degree of accuracy 
(Spearman rank correlation coefficient and p-value in the top 
left). For the HP24 null model, all frequencies are such that 
fq ^7. As such, phenotypes that are found locally are found 
only once {(j)qp = 7) and most are not found at all (the many 
downward triangular points). For the RNA 12 null model, the 
frequency of phenotypes used for phenotype p span the range 
of all fq ^ ^ (red and green) and to some phenotypes having 
fq^l (blue). As a result, we see the red and green phenotypes 
follow <fqp = fq Strongly, while the tail of the blue phenotype 
has fiuctuations between the three behaviours. 

The results from the null models demonstrate the accuracy 
of the above outlined intuition for a null relationship between 
the local connectivities of phenotypes with respect to the global 
abundance. With this in mind, we can now draw direct com¬ 
parisons between each phenotype in the null model and actual 
behaviour exhibited in the biological GP maps. For each of 
the GP maps, we plot the same phenotype as in the null model 
case. For RNA 12, positive correlations are found for each phe¬ 
notype, with deviations from (j)qp = fq being more pronounced 
for lower frequency phenotypes (blue is subject to much greater 
fiuctuations than red). The fiuctuations are approximately up to 
an order of magnitude either side of the (j)qp = fq. We see a 
similar behaviour for 82 , 8 ^ with the largest fiuctuations exhib¬ 
ited for the low frequency blue phenotype. 

Finally, we consider the biological HP24 GP map. As was 
the case in the null model version, every phenotype (apart from 
the deleterious phenotype) lies in the region where C 7. 
The notable difference in the actual GP map is the tendency for 
phenotypes with a larger fq to also be more likely to be locally 
present (log <fqp oc log fq). We can understand this with the 
following rationale: due to the neutral correlations present, if 
a single genotype with phenotype q is found locally, then it is 
also likely that other genotypes with phenotype q will be lo¬ 
cal to genotypes with phenotype p. And due to this effect be¬ 
ing more pronounced for phenotypes with a greater frequency 
iPp oc log fp, c.f. Fig. we also see this effect locally with in¬ 
creased fq resulting in a greater (fqp, leading to the positive pro¬ 
portionality between log (fqp and log fq in the actual GP maps 
when compared to the null models. 
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FIG. 10. The relationship between (pqp and fq is more complex in smaller GP maps. Local frequency of phenotypes q around genotypes with 
phenotype p (fqp) are plotted against the frequency of phenotypes q (fq) for each biological GP map and a random null model counterpart. The 
black dashed line is fqp = fq. The dotted lines are fqp = 7 and /g = 7 (c.f. Eq. ([^). In each case, three different phenotypes p with different 
frequencies in the GP map are considered (represented with red, green and blue from largest to smallest fp). (A, C and E) We plot local against 
global frequency for the random null models. 5'2,8 and HP24 illustrate the two regimes where /g > 7 and fq < 7 in all cases respectively. The 
former has local frequencies strongly determined by global frequency {fqp = /g), while in the latter, occurrences of phenotypes are rare; they 
may not occur at all (downward triangular points, <fqp = 0) or they simply occur a single time {fqp = 7). In the RNA12 null model, we see the 
blue phenotype crossing the threshold with some phenotypes having fq 7. (B, D, F) The three phenotypes are considered in each biological GP 
map. For larger frequency phenotypes (red and green in RNA12 and 82 , 3 ), we find that local frequency is, to first order, well determined by the 
global frequency in line with the random null models (up to an order of magnitude variation in local frequency in comparison to global frequency). 
For lower frequency phenotypes (blue in RNA12 and 82 , 3 ), we see that phenotype correlations are more important, an intuitive result given the 
genotypes of p will be less encompassing of the whole GP map in these cases. In HP24 all frequencies are well below the gamma threshold but 
we still see a positive (although weaker) relationship between local frequency and global frequency (unlike in the null model, where <fqp remains 
fiat with respect to fq for fq <7). This is due to the presence of neutral correlations, an effect discussed in greater detail in the main text. 








































